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ABSTRACT 


The main problem of artificial satellite theory is a restricted two body problem 
in which the Legendre Polynomial representation of the cylindrically symmetric 
potential contains only the first two terms. A generalized asymptotic expansion is 
used to obtain a first order approximation. The solution at the critical inclination 
is seen to be of a different type than at other inclinations. The solution is finite for 


all eccentricities and inclinations when suitably restricted in time. 
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NOTATION 


right ascension in spherical coordinates 

declination in spherical coordinates 

vernal equinox 

COS 19 

eccentricity of the initial instantaneous ellipse 

initial magnitude of angular momentum 

inclination 

initial inclination 

coefficient of second Legendre polynomial in the representation of the 
gravitational potential 

3/2(J2R?/ Fo) 

undetermined constants used in the solution 

h2/GM where G is the gravitational constant and M is the mass of the 
primary body 

equatorial radius of the primary body 

magnitude of the position vector to the satellite 


SIN tg 

kinetic energy of the satellite 
time 

initial value of time 

Po/ 


gravitational potential 

initial value of the argument of perigee 

angle measured from the ascending node to the satellite 
within the reference plane 

initial value of 0 


longitude of the ascending node measured from ¥ 
initial value of 2 
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A. INTRODUCTION 


The main problem of artificial satellite theory is a restricted two body problem 
in which the Legendre Polynomial representation of the cylindrically symmetric 
potential contains only the first two terms. Deceptively simple in statement, the 
main problem continues to evade satisfactory solution. We make no attempt to 
survey the vast literature pertaining to the main problem, nor to recount the history 
of the so-called critical inclination. In fact, every effort is made to make the solution 
as accessible to the non-specialist as to the astrodynamicist. 

Our goals in attempting to obtain a solution are twofold. First, we desire a 
solution employing a coordinate system based upon physical events rather than an 
abstract set of transformations. For instance, the use of averaged quantities such as 
the “mean orbital plane” or other such non-physical artifices is eschewed. Second, 
we seek a solution which does not tend to infinity at any inclination and which 
places no constraints on any dependent variable. 

Needed for the solution are a coordinate system, asssumptions concerning 
the dependent variables involved, and an algorithm for calculating the unknowns. 
These are outlined in detail in the first chapters, after which the solution is obtained 
through extensive use of MACSYMA. A discussion of the critical inclination problem 
follows. 

Because the purpose and scope of this work are nearly identical to those 
of Snider (Ref. 1], there are many similarities in organization, notation, and lan- 
guage. The principal difference in treatment lies in the geometry used rather than 


in approach. Snider’s original approach is succesfully employed using an alternate 


geometry, and his influence in this work is pervasive. 








B. ORBITAL KINEMATICS 


Figure 1 shows the reference system of spherical coordinates (r.a..3). The 
radial distance r is measured from the center of the planet O to the satellite S. The 
line Oy is in a direction fixed with respect to an inertial coordinate system. The 
right ascension a is the angle measured in the planet's equatorial plane eastward 
from the line Oy. The declination or latitude 8 is the angle measured northward 
from the equator. The position vector r of the satellite in the spherica] coordinate 
system is 


r=r(cosacos 3j}b, + r(sinacos 3)bz + r(sin 3)b3 (1) 


where (b,.b2.b3) are orthonormal base vectors fixed in the directions shown. 


polor axis 


ooo 
-* 






equatorial 
plane 


Tr 


Figure 1: Spherical coordinate system. 


Figure 2 shows how we can locate the satellite by its polar coordinates (r. 6) 
within a rotating orbital plane that contains its position and velocity vectors. Here 
8 is the argument of latitude, i.e., the angle measured in the orbital plane from the 


ascending node to the satellite. The orbital plane is inclined at an angle 7 to the 


2 











equatorial plane and intersects the equatorial plane in the line of nodes. making an 


angle 2 with the O4 line. 
orbitol plane 
“ 






equatorial plane 


Figure 2: Orbital plane 


We introduce another orthonormal set of basis vectors (B,.B;.B3) which 
move with the satellite so that B, is in the direction of the position vector r, B; is 
also in the orbital plane and B; = B, x B,. 

The basis (b,.b:.b3) may be transformed into the basis (B,.B>.B,! by a 
succession of three rotations. First the basis (b,,b;,b,) is rotated about the by 
direction by the angle . next the basis is rotated about the new first coordinate 
vector by the angle 7, and finally the basis is again rotated about the new third 
coordinate vector by the angle 6. The two sets of base vectors are related by the 


product of the rotation matrices representing each successive rotation: 


to 


cos@ sinf 0 1 0 0 cos sin 0 b; 
= ~sin6 cos@ 0) 0 cost sini ~sinQ cos? 0 bz (2) 
0 0 1 QO -sint coer 0 0 1 b3 





or 


~sin@cosQ—cos@cosisinQN —sin#@sinN+cos@cosicosQ cos@sini 
sintsin Q —sinicosQ cos i 


( cos @ cosQ — sin @ cosisinN cos @ sin + sin # cosicosQ sats: | (» 


The position vector r has only one component in the rotating basis: 
r=rB, (3) 
Using the first of Equations (2), we obtain the components of r in the fixed basis: 
r = r(cos 8 cos 2 — sin@ cosisinQ)b, 
+ r(cos@sinQ + sin @ cos? cos 2)b2 + r(sin Osin7)bs (4) 


Equating the components of Equations (1) and (4), we can obtain the following 
relations among the angles (a, 3) of the spherical coordinate system and the astro- 


nomical angles (7.1. 6): 


— 
cr 
— 


sin 3 = sin@sini 
cos 3 = cos@sec(a — 2) 


The velocity dr/dt of the satellite is obtained by differentiating (3) with 


respect to the time ?: 


dr dr dB, 
dt = We + Pay (6) 
Since the orbital plane must contain the velocity vector, we have to enforce 
dB, 
— -B;=0 7 
qi 3 (v) 


Substitution of Equation (2) into Equation (7) leads to a relationship which uncou- 
ples the equations for (6) and 7(@): 
dQ tan di 


—— 


dO sini dO 


io 4) 
~~ 


4 





The velocity (6) can then be written 


dr dr dé di 
ae Press + ie ( + tan @ cot 5) B, (9) 


In the following part of this paper, we obtain expressions for r(@), 2(@), 2(8), 
and dt/d6(@). The position and velocity vectors of the satellite are then calculated 
from the formulas in this chapter. The classical orbital elements p,e, and w are 
the semilatus rectum, eccentricity, and argument of perigee of the instantaneous 


{osculating) ellipse determined by the position and velocity vectors. If needed, p(@), 


e(6), and w(@) can be obtained from our solution r(6) and dt/d0(@): 


Numerical integration of the expression for dt/d0(@) allows a transformation between 
6 and t. Figure 3 shows the satellite at time ¢ and its relationship to the various 
variables used. 

This method employs the true, rather than mean, orbital plane to specify the 


satellite’s position and is due to Struble (Ref. 2, 3, 4]. 


C. EQUATIONS OF MOTION 


The expressions in spherical coordinates for the kinetic and potential energies 


per unit mass of a satellite orbiting around an oblate planet are respectively: 


1 | fdr\’ dB\? da\? 
r=3{($) +r? (F) + teoe'a() | (10) 








; 2 
ye oe t+ 5 (1 asin") (11) 
5 





Figure 3: Satellite at time 1. : 


where G is the gravitational constant, Af is the mass of the planet, and R is the 
equatorial radius of the planet. Substitution of these equations into Lagrange’s 


equations 


ddaT-V) 0 nL = 
di a(a@y a7) = g=7,a, or B 


results in the following equations of motion: 


dr d3\* »>,(da\? a 
= sons = ) 

dt? (Fr) Poe (3) Or (12) 

df. 2,da\ _ 

a (+ cos oF = 0 
d { .d3 2 ,{da an OV 
77 (> at? sin 3eos 3 (Ft) ~~ O93 (13) 

6 





Initial conditions are established by requiring that at the initial time to the 


orbital parameters of the Keplerian two body ellipse, determined by the initial posi- 
tion and velocity vectors, are known. The actual orbit is then tangent to this initial 
instantaneous ellipse (see Figure 3). Equating the initial position and velocity vec- 


tors given by Equations (3) and (9) to the two body expressions, we obtain 


Oo) = Tea Fle) = Say) 
F (to) = ee ee (15) 
r3 [1 + tan 6 cot in $(G0)| 
0) =% (16) 
(89) = NA (17) 


Here ho = /GMpo is the initial value of the satellite’s specific angular momentum 
about the center of the planet. The subscript 0 on a symbol denotes that the 
parameter is evaluated at the initial time fo. 


We immediately have two integrals of the equations of motion: 
T+ V = constant (18) 


and 


d 
2S constant (19) 


2 

r’ cos’ B— 

Equation (18) states that the mechanical energy of the satellite remains constant. 
Now, from Equations (1) and (19) 


dr 2 


d. 
7 Ma a cos? I = ho COs ip (20) 








Equation (20) states that component along the polar axis of the specific angular 
momentum of the satellite remains constant. Inserting Equations (3) and (9) into 


Equation (20), and applying the initial condition (15), we obtain 


dt = r* cosi di 
7 cous (1 + tand ots) (21) 


This allows the independent variable to be changed from t to @. 





Letting u = po/r, and using Equations (5), (20), and (21), we can rewrite the 


remaining equations of motion (12)-(13): 





di —2Jusin 6 cos @ sini cos* i (22) 
do =< + 2Jusin? 6 cos? i 

du cos*2 Jcos*i[ 4 scp ae 

as a2 [u (1 — 3sin* @sin*z) 


+ 2uM" sin 0 cos6(1 — Scos?s) — dua" sin? @cos?i-2(%) sine | 23) 
a6 sin U cos cos’ 2 ume sin cos? 1 a Sl cos” 2 ( 


2, eine 6 
4J*u sin” 8 cos a ~ cos (2+ sin’ i) + 5 (u a) sin 6 cos* ] 


The terms with d?u/d6? on the right side of (23) can be eliminated, yielding the 


ci 


equivalent equation 


d*u 


do dé 


+u= {cPcosti ~ Jet u COs *i [25 sin 8 cos (3 cos w—1)- 3 
_ wre sin? 8 cos 8 cos® i(3 — cos? i) (24) 


/(c4 + 4Juc? sin? 6 cost i + 4u?J? sin? 6 cos’ i) 


Here we have introduced the shorthand notation c = cost, § = Ssinto, J = 








The differential Equations (22)-(23) are coupled by the nonlinear terms and 


apparently cannot be solved analytically. If we expand the right sides of (22) and 


(24) in a Taylor series expansion in powers of J, the equations simplify to 


di —2Jusin@cos @sinicos*i 


+ 4J?u?sc’ sin? 6 cos 8 + O(J%) (25) 





do (od 
du cos? Jcos?i(~—4usin?@costi = ie 
7 tue a2 4 OE WL + sin 6(7 cos” i — 3)] 
du . 2: du - 2 2. 9 
+ 2u sin 6 cos 6(1 — 3cos*i) — 2 (3) sin* 6 cos iS (26) 
2 «2 6; : 2 4° 
4J = Lut + 3sin? 61 — 2cos?i)] + 3u sin a 2 
c c 
du . . 2 . - 2 2 . du si . 2 2 . 3 
+u; sin 8 cos 6[4 sin 2+ sin* 0(1 — 3cos”2)] + 9) sn 6 cos iS + O(J*) 
Using (8) together with (25) yields 
dQ +2 9 3 a4 3 
70 = —2Jucsin® @+ 4J*u’c’ sin* 8 + O( J”) (27) 


Each of the neglected terms in Equations (25)-(27) are indicated by the O 
symbols. These are terms which will be multiplied by J to the third power or higher. 
Note that the Equations (25)-(26) are identical to those used as the starting point 


in the analysis of Eckstein, Shi, and Kevorkian [Ref. 5]. 


D. METHOD OF SOLUTION 

We will use a perturbation technique to solve Equations (25)-(27). Following 
Erdélyi [Ref. 6], we define the order relations O and o as follows. For two real-valued 
functions f(r) and g(x), we define the relation f = O(g) if there exists some real 
constant € such that |f] < |g] for all z in some domain of interest. Similarly, we 
define f = o(g) as r — 7, if for all € > 0, there exists some neighborhood of x, such 


that |f| < elg| within that neighborhood. 


9 











The series >, a,(0, J) J” is said to be a generalized asymptotic approximation 


to N terms of f(@,J) with respect to the scale {J"} as J - 0 if 
N 
f(0,J) = ¥ a,(0, J)J™ + o(J%) as J 0 
n=0 


Let D be some subset of the real line to be determined. We will say that the 

above generalized asymptotic approximation holds uniformly for 6 € D if 

N 

f(6,J)->> a,(0,J)J% = o(J®) as J 0 

n=0 
uniformly for 6 € D. In order to get uniformity it will sometimes be necessary to 
bound the independent variable @, thus determining D. Note that these definitions 
differ from the usual asymptotic expansion definitions in that we allow the coefficient 
functions a; to be functions of both @ and J. 


Let us suppose that each of 2, u, and 2 have formal generalized asymptotic 


expansions for some suitably restricted values of 0: 


it 3 i,(0, J) J* (28) 
k=0 

Une = u,(9, J)J* (29) 
k=0 

Qs > 0,(6, J) J* (30) 
=0 


Brenner and Latta [Ref. 7] found that a modification was needed to their 
time-like variable M in order to avoid resonance. They accomplished this by means 
of an additional variable w which multiplied M, allowing secular term elimination. 
Analogously we introduce the variable y which may be added to integer multiples of 


6. In this way we too will avoid resonances except perhaps at certain inclinations. 


10 








The explicit form is then 


y= 3 (6, J)J* (31) 


k=0 


We further stipulate that for J = 0, u assumes the Keplerian two body 


solution and that y is the true anomaly. This forces 


y= O-w, +Jy+J*y24+--- (32) 
true anomaly 

Un 1+ eocosy +Ju, + J*u,4--- (33) 
a / 


Keplerian solution 


An algorithm for the perturbation procedure is then: 


Let n=1 

Substitute Equations (28), (30), (32), (33) into the equations of motion 
Equate coefficients of J” 

Solve for the n“* order solution 

Iterate on n 


Having chosen a coordinate system, made assumptions concerning the vari- 
ables, and given an algorithm for determining the unknowns, we are prepared to 
give the solution in the following chapter. 

E. SOLUTION 

Substituting Equations (33) and (28) into (25) and equating the terms mul- 

tiplied by J yields 


diy 


0 7 —scsin 26 — ee sin(y + 26) + = sin(y — 26) (34) 


A solution to this equation is 
11 = S cos 26 + “F cos(y + 20) + FS cos(y — 26) 


+ Kcos(2y — 26) + K, (35) 








The last two terms may be added because they are, to lowest order, homogenous 
solutions to Equation (34). The term multiplied by the constant A, was added to 
eliminate secular terms in 22. Note that differentiating this term with respect to @ 
produces terms multiplied by J, from Equation (32). The constant K, was added 


to satisfy the initial condition (16), which implies that 7;(4) = 0 so 


Kk, = -F cos 26 — see .95(3 — Wo) — == cos(9o + wo) — Ay cos 2wo 


6 
Substituting Equations (28), (32), (33), and (35) into (26), and equating 


terms multiplied by J yields 








du; 3s? 2 5s? 1 2\ 2 2 
Ge test- Sg (S 41 + qI(2 + 5¢9)s — 2e6] cos 20 

€ 0 2 

+ 7 (-9s? + 8) cos 2y + 3 (lls — 6) cos(y + 20) (36) 
2 2 , 
+5236? — 2) cos(2y + 20) + [Bes ap) 2 a cos(2y ~ 20) 
2sK2 dy, : d?y, , 
: eo (2h +4-5s cosy + €o—75 sin y¥ 


In the above equation, the cos y and sin y terms would produce secular terms @ sin y 
and @cosy in u,. Choosing dy,/d? = 5s?/2 — 2 will eliminate these possibilities. 
Integrating yields 


2 
y= (= = 2) (6 = 9) + K[sin(2y ~ 20) + sin 2wo| (37) 


The term multiplied by 3 was added to eliminate secular terms in uz. The constant 
terms in (37) were added to satisfy the initial condition y(@)) = 09 — wo. 
A solution to lowest order of Equation (36) is then 


2 8 oe 1 
uw=l- = +e (== + 7 os —s?(2 + 5e5) + 2€2] cos 26 


2 
+33 (98" — 8)cos2y + sp(c ls? + 6) cos(y + 26) 
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2 2 
+ 5o(—3s? + 2) cos(2y + 28) + $38? — 2) 


- 2sky 
c 


cos(2y— 20) (38) 


_ esky 
c 





+ K4cos(y — 26) 
+s cos(y — 99 + wo) + Kesin(y — 09 + wo) 


The term multiplied by K4 was added to eliminate secular terms in uz. The terms 
multiplied by AK, and Kg were added to satisfy the initial conditions (14). 

The calculations proceed by substituting Equations (28), (32), (33), (35), 
(37), and (38) into (25) and equating terms multiplied by J?: 


diz sce2(15s? —J4)] . 
————_. - 2y —26)+4+... 
do (ped) [ee (38) 


For brevity we have indicated on the right side of Equation (39) only the term that 
would produce secular terms in 72. Removal of this term by making its coefficient 
zero determines Ay. Equation (39) is then integrated to determine 2. 

It should be noted that 22 is needed for the determination of the constants 
K3 and K4 in the equation for ug, so integration is required. 

Continuing the procedure by equating the terms multiplied by J? in the 
expansion of Equation (26) determines y2, A3, and Ay. Final values of all the 
constants are listed in the Appendix. 

00(8) is determined by substituting (30) into (27) and proceeding in the same 
way as above. Note that terms in J?6 arise in (0). These must be retained and 
will restrict our variable 6 accordingly. 

In the form below, use has been made of trigonometric identities in order to 
isolate terms containing the quantity (5s? — 4) in the denominator. It may be seen 
that when (5s? —4) is zero, each of the variables has a finite limit. Note the necessity 


of keeping the J70 terms in y and 9. 











When the quantity (5s? — 4) is zero, the terms multiplied by J and J?6 


combine, allowing the division by (5s? ~— 4) to occur. 


5s 


r= po/{1 + eocosy + J[1— 2S +63 (i- =) 


+ 5(- + 5e2)s? + 2e?) cos 20 + (93 — 8) cos 2y 


+ s(- 11s? + 6) cos(y + 26) + (_ 3s? + 2) cos(2y + 26) 


2 
+ 238? — 2) cos(2y — 26) 


eo[15(2 + e2)s4 — 14(4 + €2)s? + 24] sin [70 (= - 2)| sin(# + wo} 





i 12(5s? — 4) 
e2s?(15s? — 14) sin [Jo es - )] sin [2w0 —J6 (= ~ 2)| 

- 6(5s? ~ 
eps” € 9.2 

= cs cos(y = % + 3u9) + 94 (38 = 2) cos(y — 365 + dwg ) 

e2s? 

- < cos(y — 509 + 3wo) + (38? — 2) cos(y — 28 + 2wo) 
3€0S : €o, 2 

some cos(y — 449 + 2wo) — 7s +1) cos(y + 2wo) (40) 
if 

+ gl(-2 + 5e2)s? — 2e?) cos(y + 0 + wo) 


] 

+ 7[6 + 5¢6)s’ ~ 4(1 + €6)] cos(y — 80 + wo) 
1 

+ 3 + 5e3)s? + 2€?] cos(y — 309 + wo) 


2 
+ ia — 4) cos(y + 389 — wo) + P(-7s? + 6) cosly + 09 — wo) 
+ $8 (54? + 4) cos(y ~~ wo) + + (2s? — 1) cos(y + 29) 
+ + (-38? +1)cos(y — 269) + + (-3s? +2) cos y 

2 €9s? 
+ egs” cos(Op + wo) + a cos(36> — wo) 


+ 8 cos 266] n O(J?, 1°6,...)} 





14 


4— on + J{ (= -2) (8 — 6) 


e3(—75s° + 260s* — 296s? + 112) sin [J6 (8% — 2)| cos [2wy — Jo (5 — 2)| 


24(5s? — 4)? 
2 <2 2 2 
anf €487(—15s? + 14)(15s? — 13) cos2wp  egs? 
J 0 ey a (5s? — 13) cos(6 + wo) 
C08 jo 9 ree 
—G (15s — 13) cos(309 — wo) + 5 (15s — 13) cos 269 (41) 


sgl 5(9e2 + 34)s4 + 4(Ge2 — 34)s? - sed) + O(J?, J°6,...) 


= ot sed {5 cos 26 + > © cos(y + 20) + — = cos(y — 20) 


e2(—15s? + 14) sin [Je (= - )] sin m — J6 (3 _- 2)] 
12(5s? — 4) (42) 


~ 5005 20 - F £08(3% —w9) —- 5 £08( 4 + wo) + O(J?, J76,...) 


Q = Qe + 0 {Bo ag x sin 20 aay 


+ = sin(y +20) ~ 5 sinty — 26) — sin 2 


4 egain( Oy = tag) = S sin( 30 = 5 sin( + wo) 
e2(15s4 — 45s? + 28) sin [ve (= - )] cos [2w0 — J60 (32 —_ 2)| 
— O6(HSR=AYEO—~—<C—sSOSSC‘CS 


oJ*0 egs"(15s" — 14) cos 2w9 — €98” cos(o + wo) (43) 


12(5s? — 4) 


re 
- < cos( 34 — wo) ~ s? cos 28 + % (7s? — 4) 


+ aI —s *+6)} + 0(J7,J%,...) 
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@ ~3527 49 
t= t+ f ft 4 s[ PEF?) cos a0 
9 


ho 
—45? 4+ 3 
+ €9(s? — 1) cosy + cot costy + 26) 
9,2 
4 col tO costy — 20) (44) 
e3s?(15s? — 14) sin [J6 (§° — 2)] sin [29 — Jo (8 - 2)] 
is 12(5s? — 4) 


2 2 2 
fe opt > cos 26 + a cos(38) — wo) + — cos( 69 + 0)| 


4 O(J?, J%0,...)}d0 


To check the solution, we can see if the specific mechanical energy (18) of 
the satellite remains constant. Substitution of the solution (40)-(42) into Equations 
(10) and (11) yields 


GM(1—¢€3) _ GMJ,R?(1 — 3sin? 3p) 


T+V=- 
2Po 2[r(to)} 


HOS I 


The first two terms on the right side are easily recognized as the value of the specific 
mechanical energy at the initial time fo. 

As a further check on the solution, we can see if it reduces to previous results 
for equatorial] and polar orbits, obtained by Danielson and Snider [Ref. 8}. Setting 
io = 0 and using the independent variable a measured from the line O04, we find that 
Equations (40)-(42) reduce to Equations (18)-(22) in [Ref. 8]. Setting 79 = 7/2 and 
using the expansion cos(y+Jk) = cosy—Jksiny, we find that Equations (40)-(42) 
reduce to Equations (38)-(41) in (Ref. 8].' 


1There are misprints in Equations (37)-(40) of [8]. The term 2eqc2 cos y should be added on the 
right of Equation (37). The third from last term in Equation (38) should be e9/3 cos(3.39—wa). The 
term —c2J?@ should be added on the right of Equation (39). The sign of the second trigonometric 
term in the expression for cg should be changed: (—e9/12 + €3/32) sin( do + wo). 








F. THE CRITICAL INCLINATION 


There exist many “solutions” to the main problem. Most of these are un- 
satisfactory. Some are so abstruse as to have little or no physical meaning, while 
others suffer from having a restricted domain in inclination, eccentricity, or both. 
Chief among the difficulties with proposed solutions is their behavior at certain in- 
clinations. The most persistent of these troublesome inclinations is the so-called 
critical inclination toy, = tarcsin2//5 (mod mr). This manifests itself in the form 
of denominator terms like (4 — 5sin?z), (1 — 5cos?z), or (tani — 2). When 7 as- 
sumes this critical inclination, zi, each of these terms become zero, rendering the 
solutions useless. It is at this point that we must confront the divisor (5s? — 4) in 
Equations (40)-(44). As remarked, each of (40)-(44) has a finite limit at the critical 
inclination. How does this come about? At the critical inclination, the quotient 
sin( J0(5s?/? — 2)/(5s* — 4) is replaced by the limit, yielding a term in J@. This, in 
the parlance of orbital mechanics, is an odious secular term — an unbounded term 
which grows with the time-like variable. The expressions (40)-(44) at the critical 


inclination are: 


3s? 5s? 
C= pol{1 + eacosy + Jfi~ 9 4g (1-*) 
] ‘ re 2)\.2 ) 9,2 : eG 2 «9 : 
+ + 75(-© + 5e)s + 2e,) cos 20 + hag — 8) cos2y 
€ 2 € 2 
a og tts + 6) cos(y + 20) + 5 )(—3s + 2)cos(2y + 28) 


2 
+ F (3s? — 2) cos(2y — 26) 


ezs? ms) 

- aT cos(y ~ 09 + 3wo) + risa — 2) cos(y — 349 + 3u9) 
e2s? €0 

Pt ares cos(y — 509 + 349) + 7 (38? — 2) cos(y — 269 + 2uo) 


3€98? 





€ 
cos(y — 46y + Quy) — ais? + 1) cos(y + 2u0) (45) 
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1 
ea gl(-2 + 5e2)s? — 2€2] cos(y + 00 + wo) 


1 
+ ql(6 + 5e2)s? — 4(1 + €3)] cos(y — 9 + wo) 


1 
+ Zl-(d + 563)s? + 242] cosly — 300 +o) 


2 2 
+ 33 (98? — 4) cos(y + 30 — wo) + q(-78? + 6)cos(y + % — wo) 
a 
+ ae —5s? + 4) cos(y — 0) — me = (28? — 1)cos(y + 24) 
4 f(—3s? +1) cos(y — 20) + 24 +2) cosy 


22 
+ €98* cos(O + wo) + cos(38 — wo) + s* cos 206] 


+ ols (15(2 + €2)s — 14(4 + €2)s? + 24) sin(@ + wo) 


+ SO (150 - 14) sin(2u9)| + O(J?, J°6,.. y 


h a) 
30 { Fq(—105s4 + 130s? — 28) cos wo + (153? — 13) cos(4 + 0) 
€98? 5 $2 
gids — 13) cos(30@ — wo) + > 


a [5195 + 34)s4 + 4(9e2 — 34)s? — sed] | + O(J?, J76,...) 


(15s? — 13) cos 20) (46) 


1 io+ sed {5 cos 20 + = cos(y + 26) + > cosly — 28) 


- 5 £05 26 — cos( 3d ~ we) — 2 cos (Mo + wo) } (47) 
+ prose 15s? + 14) sin Quy + O(J?, J*6,....) 
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Q = M% + ed {Bo — 8 + 5 sin 20 — ensiny 
€o . € . oe 
+ sin(y + 26) — 2 sin(y — 20) — 5 sin 265 
+ egsin(@ — wo) — = sin(36 — wo) — > sin(60 - wo) 
2 
+ oJ*0{ 2,65? — 7) cos(2w») ~ €gs? cos( + wo) (48) 
€ 


2 
€9S 
- ao cos(309 — wo) — s? cos 209 + ritka — 4) 


mn a(-3 + 6)} + O(J?, J6,...) 


6 —357 442 
t = t+ | fi ¢ [SPS *9) cos 26 
0 


ho 
“set 
+ €9(s? — 1) cosy + cot costy + 20) 
25? 
sp SOE Ee eg (49) 


9 


5 wees 4 Eos? Eps? 
+ s°—-1l4+ "5 COs 269 + ar cos(389 — wo) + a cos(99 + wo) 


242 
POL (15s? — 14) sin(2w9) + O(J?, J°6,.. jaa 


+ 


As can be seen, none of the variables fails to be defined at any inclination or 
eccentricity. 

A long-standing debate in astrodynamics centers about the nature of the 
critical inclination’s ubiquitous presence in proffered solutions. Of course something 
important does occur at the critical inclination — the line of apsides remains fixed. 
That this event should also cause most attempts at solving the problem to fail has 
been the subject of much debate and, apparently, misunderstanding. 

Examination of Equation (45) shows that when 7 assumes the critical incli- 


nation there result terms in J?0. At the order to which we have approximated the 
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solution the form of our approximation at the critical inclination differs from that of 
all inclinations. In fact, if 7 # tj, it appears that we can obtain arbitrary accuracy 
for sufficiently small J and @ in some J dependent interval [0,c/J™}, (m constant, 
integer, c an O(1) constant), by taking sufficiently many terms. However, when 
t = tqit, the expansion has a restricted range of validity, for 6 in [0,c/J] where c is 
an O(1) constant. In other words, while we have dropped the higher order terms and 
higher order secular terms in our solution, it is hoped that away from the critical 
inclination it will be possible to improve the approximation to arbitrary accuracy 
simply by continuing the method used. Not so at the critical inclination—we are al- 
ready confronted with a secular term which cannot be eliminated, regardless of how 
many terms we keep in our approximation. We are thus in agreement with Coffey, 
Deprit, and Miller (Ref. 9] that the critical inclination is an intrinsic singularity, 
independent of the method used to solve the problem. Moreover, there may emerge 
other critical inclinations as higher order analyses are performed. It may be that 
a higher order approximation will yield other inclinations where secular terms arise 
which cannot be eliminated. The well-known critical inclination might be only the 
first in a sequence. 

Enters now the astrodynamicist, the man in the field. When the equations 
of motion are integrated at the critical inclination, no strange effects are reported. 
The mathematicians have declared the critical inclination a glaring singularity, yet 
everyday practitioners consider it a trifling locale. If a satellite were to start an 
orbit at the critical inclination, Equation (45) predicts a secular variation in u(6) 
of order J?6. Clearly such slow growth with @ explains why no appreciable effect of 
the critical inclination has been encountered in the intensive numerical integrations 
which have been carried out over the years. Additionally, if a satellite is initially at 


the critical inclination, Equation (47) states that the expression for its inclination 
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also contains a term in J*0, This means that the satellite will tend away from the 


critical before the secular terms which arise in u(@) at the critical become important. 
As the satellite moves away from the critica] inclination, the terms in u(#) resume 
their bounded trigonometric form. From a practical viewpoint, the critical inclina- 
tion problem is in fact a nonproblem. The simplifications and assumptions made 
will perforce limit the time over which the solution is valid, further obscuring at- 
tempts to physically observe anomalies at the critical inclination. We must concede 
that other perturbative forces would play a major role in determining the degree to 
which the critical inclination could induce observable effects. This is not to endorse 
a variety of patchwork measures now being used to determine satellite motion, such 
as “gapping” the critical inclination in order to avoid dividing by zero. The critical 
inclination, delicate as though it may be, is perhaps a powerful telltale for faulty 


modeling. 


G. CONCLUSIONS 

As shown, our solution satisfies the goals given in the introduction. 

If the perturbation procedure is continued, it is anticipated that each of the 
additional terms added to Equations (40)~(44) will be multiplied by one of the 
factors (J*, J*, J30, J40,...). If we restrict 6 < 1/J, the neglected terms should be 
order J?. (For an Earth satellite J < 3/2 x 107%, so for at least 100 revolutions the 
relative error should be order 107°.) 

If we restrict 8 < 1, all of the terms of order J?0 in the solution (40)-(44) can 
be dropped without increasing the order of magnitude of the error. This considerably 


simplifies the solution to: 
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5s? 
po/ {1 + egcos |9 — wo + J(8 — 4) ge 2 
3s? ae 1 ie nie 
t - ge 1- — + 75 — (2 + 5ea)s + 2e5) cos 20 
$8 (942 — 8) cos(26 — 2w9) + sa(- 11s? + 6) cos(30 — wo) 
di + 2) cos(46 — 2wo) + £2,362 — 2) cos 2wo 
e2s? e 
a “Tg 0089 - % + 2Qwo) + 5g 8S — 2) cos(@ - 39 + Quo) 


i 
7 coe 505 + 2wo) + (3s? — 2) cos(O — 249 + wo) 


3e9s7 





cos(@ — 44 + wo) — Fs? + 1) cos(@ + wo) 
al(-2 + 5e2)s* — 2e?] cos( + 9) 


AG + 5e2)s? ~— 4(1 + €2)] cos(@ — 8) 


| 


om —(14 + 5e2)s? + 2e2] cos(@ — 34) 


(99 24) cos( 8 305 = 006) + $2(_768 + 6) cos(8 + 0 — 2w) 
(59 + 4) cos(@ — 8 — 2wWo) + + (2s? — 1) cos(@ + 269 — wo) 
+ (-3s? +1) cos(@ — 209 — wo) 

+ (-3s? + 2) cos(9 — wo) + €08” cos(4o + wo) 


€98” 


qa 008(3%0 — wo) + 8? cos 266] + O(J?, J°6,.. )} 


io + sed 5 cos 26 + > © cos(30 — wo) + 5 c08(4 + wo) 


5005 209 — F£08(300 — wo) — 5 00s(4o + 0) + O(J?, J36,...) 


to 
bo 


Qo + cJ 6 —O+ sin 20 — egsin(8 — wo) + = sin(30 — wo) 


1 
> sin(@ + wo) — 5 sin 209 + eg sin(99 — wo) — > sin(309 — uo) 


2 sin(Oo + w0)| +0(J?, J°0,...) 


8 9,2 
= t+ {14 J[ PSF? cos 
ho Cr 2 


€9(—4s? + 3) 


+ €o(s? — 1) cos(@ — wo) + cos(30 — wo) 


2 
cos(9 + wo) + 3? —14 > cos 2 


€9(—2s? + 1) 
- 2 
€9s? €oS? 
+ 7 cos(34@ — wo) + 5 cos(6> + u»)| 


‘ o(s?, 3°, ...) bao 
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APPENDIX 


__ eseg(—15s? + 14) 


K,= 
: 24(5s? — 4) 
é 
ky = -> cos 20 — = cos(30y — wo) — > cos(O + wo) 
2 2 
cse3(15s? — 14) 
Shol —  cos tu 
Misst=4) 
Hex e2(—75s° + 260s* — 296s? + 112) 
oa 48(5s? — 4)? 
K,= og tile + 2)s4 — 14(e2 + 4)s? + 24] 


5s? — 4 


2 2 
K; = =3(—98? + 8) cos(20) — 2w) + 5 (38? — 2) cos(49y — 2wo) 
~ (es? + 4) cos(9 + wo) + (3 — 2) cos(309 — wo) 
2 
+ F(-3s? + 2) cos Quo — 5602 — €2)s? + 2e2] cos 24 


1 
+ 79 15¢0 + 18)s? ~ (e3 +1) 


2 2 
Kg = +2(6s? — 5)sin(209 — 2wo) + <3 (387 + 1)sin(48) — 2wo) 
1 ; . 
+ s[eo(—s? + 1) + 2K4]sin(Bo + wo) + (38? — 2) sin(8 — wo) 
€o 2 : e2 2 . 
+ gic’s + 2) sin(309 — wo) + i ees + 1)sin 2wo 


1 
+ gl-(5e0 + 2)s? + 23] 
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[5] 


[6] 
[7] 


[8] 


[9] 
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